% To analyze the time series related to the
% Tsunami 2010 Chile

%load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\IPC_FEB2010;
%load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\PPT_FEB2010 ppt_data

load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\ipc20100301vsec.sec.mat;
ts1 = data.DATA(:,3)-nanmean(data.DATA(:,3));
fday = data.FDAY;

if sum(isnan(ts1)) > 0,
    display('NaNs');
     pause;
end;
    
%ts1(isnan(ts1)) = -3.221; (feb 17)

  b1 = min(fday):(1/24):max(fday);
  sp=spline(b1,ts1(:)'/spline(b1,eye(length(b1)),fday(:)'));
  ts=ts1-ppval(fday,sp);
    
% ts=ts1;  
%waven  = 'morl';
waven = 'cgau4';


subplot(211);
plot(fday,ts);
datetick;
axis([-inf inf -inf inf]);
perdiod_in_minutes = [1:100];

% Sc = scal2frq(1./(fliplr(perdiod_in_minutes)*60),waven,1); %to get desired scales for a set of frequencies
% Sc = scal2frq(logspace(-3,-1,10),waven,1); %to get desired scales for a set of frequencies
Sc = [100:10:1000];

c = cwt(ts,fliplr(Sc),waven);%Complex wavelet transform
perd = 1./scal2frq(Sc,waven,1);
subplot(212);

%imagesc(fday,fliplr(perd)/60,abs(c));
imagesc(fday,fliplr(perd)/60,abs(c(:,1:10:end)));
datetick;
caxis([0,12]);
colorbar('horiz');

dir_path = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\';
S = dir([dir_path 'ppt*.sec']);

for i = 1:length(S),
     display([dir_path S(i).name]);
     data = obs_read_intermagnet_second([dir_path S(i).name],'XYZF');
     eval(['save ' dir_path S(i).name '.mat ' 'data']);
end;

dir_path = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\';
S = dir([dir_path 'ipc*.sec']);

for i = 1:length(S),
     display([dir_path S(i).name]);
     data = obs_read_intermagnet_second([dir_path S(i).name],'XYZF');
     eval(['save ' dir_path S(i).name '.mat ' 'data']);
end;

dir_path = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\';
S = dir([dir_path '*.mat']);

ipc_fday = zeros([1,12*86400]);
ipc_z = zeros([1,12*86400]);
st = 1;
en = 86400;

for i = 2:length(S),
    
    eval(['load ' dir_path S(i).name]);
    ts = data.DATA(:,3)-nanmean(data.DATA(:,3));
    fday = data.FDAY;
    
    ipc_fday(st:en) = fday;
    ipc_z(st:en) = ts;

    st = en+1;
    en = en + 86400;
end;

load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_z_20100218_20100301 ipc_fday ipc_z;

ts = ipc_z(1:10:end);
fday = ipc_fday(1:10:end);

if sum(isnan(ts)) > 0,
    display('NaNs');
     pause;
end;
    
%ts1(isnan(ts1)) = -3.221; (feb 17)

  b1 = min(fday):(2/24):max(fday);
  sp=spline(b1,ts(:)'/spline(b1,eye(length(b1)),fday(:)'));
  ts1=ts-ppval(fday,sp);
    
% ts=ts1;  
%waven  = 'morl';
waven = 'cgau4';


subplot(211);
plot(fday,ts1);
datetick;
axis([-inf inf -inf inf]);

% Sc = scal2frq(1./(fliplr(perdiod_in_minutes)*60),waven,1); %to get desired scales for a set of frequencies
%Sc = scal2frq(logspace(-3,-1,50),waven,10); %to get desired scales for a set of frequencies
Sc = [10:2:200];

c = cwt(ts1,fliplr(Sc),waven);%Complex wavelet transform
perd = 1./scal2frq(Sc,waven,10);
subplot(212);

%imagesc(fday,fliplr(perd)/60,abs(c));
imagesc(fday,fliplr(perd)/60,abs(c(:,1:10:end)));
datetick;
caxis([0,12]);
colorbar('horiz');

%% cross spectra
%cross wavelet spectra
Sc = 100:10:1000;
waven = 'cgau4';
perd = 1./scal2frq(Sc,waven,1);
ipc_dir = dir('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\*.mat');
ppt_dir = dir('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\*.mat');


for i = 2:13,
    
eval(['load ' 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\' ppt_dir(i).name]);
ppt_z = data.DATA(:,3);
fday = data.FDAY;
ppt_h = data.DATA(:,1);
eval(['load ' 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\' ipc_dir(i).name]);
ipc_z = data.DATA(:,3);
ipc_h = data.DATA(:,1);
display(['D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\' ppt_dir(i).name]);
b1 = min(fday):(1/24):max(fday);

sp=spline(b1,ipc_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_z_f = ipc_z-ppval(fday,sp);

sp=spline(b1,ppt_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ppt_z_f = ppt_z-ppval(fday,sp);

sp=spline(b1,ipc_h(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_h_f = ipc_h-ppval(fday,sp);
sp=spline(b1,ppt_h(:)'/spline(b1,eye(length(b1)),fday(:)'));
ppt_h_f = ppt_h-ppval(fday,sp);

subplot(211);
plot(fday,ipc_z);
datetick('x','keeplimits');
a = axis;
title([ipc_dir(i).name(1:11) 'Z weighted with PPT X spectra']);

clear ipc_z ipc_h ppt_z ppt_h;

ipczw = cwt(ipc_z_f,Sc,waven);%Complex wavelet transform
ipchw = cwt(ipc_h_f,Sc,waven);%Complex wavelet transform
ppthw = cwt(ppt_h_f,Sc,waven);%Complex wavelet transform
pptzw = cwt(ppt_z_f,Sc,waven);%Complex wavelet transfor

Hxy = ipchw .* conj(ppthw);
cc = abs(max(max(abs(Hxy))) - abs(Hxy) );
weight = (cc./max(max(cc)));
dd = (abs(ipczw).*weight);


subplot(212);
imagesc(fday,perd/60,abs(ipczw).*weight);
caxis([0,12]);
datetick('x','keeplimits');

saveas(gcf,['D:\Manoj\tsunami\2010_Chile_Easter_Island_data\'  ipc_dir(i).name(1:11)],'fig');

eval(['save ' 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\' 'Wav' ipc_dir(i).name(4:11) '.mat ' 'ipchw ipczw ppthw pptzw'  ]);

end;

%% Cross spectra PPT

Sc = 100:10:1000;
waven = 'cgau4';
perd = 1./scal2frq(Sc,waven,1);
ipc_dir = dir('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\*.mat');
ppt_dir = dir('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\*.mat');


for i = 2:13,
    
eval(['load ' 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\' ppt_dir(i).name]);
ppt_z = data.DATA(:,3);
fday = data.FDAY;
ppt_h = data.DATA(:,1);
eval(['load ' 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\' ipc_dir(i).name]);
ipc_z = data.DATA(:,3);
ipc_h = data.DATA(:,1);
display(['D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\' ppt_dir(i).name]);
b1 = min(fday):(1/24):max(fday);

sp=spline(b1,ipc_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_z_f = ipc_z-ppval(fday,sp);

sp=spline(b1,ppt_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ppt_z_f = ppt_z-ppval(fday,sp);

sp=spline(b1,ipc_h(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_h_f = ipc_h-ppval(fday,sp);
sp=spline(b1,ppt_h(:)'/spline(b1,eye(length(b1)),fday(:)'));
ppt_h_f = ppt_h-ppval(fday,sp);

subplot(211);
plot(fday,ppt_z_f);
datetick('x','keeplimits');
a = axis;
title([ppt_dir(i).name(1:11) 'Z weighted with IPC X spectra']);
set(gca,'FontSize',16);
clear ipc_z ipc_h ppt_z ppt_h;

eval(['load ' 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\' 'Wav' ipc_dir(i).name(4:11) '.mat ' 'ipchw ipczw ppthw pptzw'  ]);

Hxy = ipchw .* conj(ppthw);
cc = abs(max(max(abs(Hxy))) - abs(Hxy) );
weight = (cc./max(max(cc)));
dd = (abs(ipczw).*weight);


subplot(212);
imagesc(fday,perd/60,abs(pptzw).*weight);
set(gca,'FontSize',16);
caxis([0,12]);
datetick('x','keeplimits');

saveas(gcf,['D:\Manoj\tsunami\2010_Chile_Easter_Island_data\'  ppt_dir(i).name(1:11)],'fig');


end;

L =ipc_fday>points3(2).Position(1) & ipc_fday<points3(1).Position(1);
selected_z = ipc_z(L);
selected_fday = ipc_fday(L);

b1 = min(selected_fday):500/(3600*24):max(selected_fday);
sp=spline(b1,selected_z/spline(b1,eye(length(b1)),selected_fday'));
spline_fit = ppval(selected_fday,sp);
L = selected_fday > points4(2).Position(1) & selected_fday < points4(1).Position(1);
selected_z(L) = spline_fit(L);
L = ipc_fday > points4(2).Position(1) & ipc_fday < points4(1).Position(1);
LL = selected_fday > points4(2).Position(1) & selected_fday < points4(1).Position(1);
ipc_z(L) = spline_fit(LL);

save d:\manoj\tsunami\curser_info_correct_ipc_data_1 cursor_info points2;
save d:\manoj\tsunami\curser_info_correct_ipc_data_2 points3 points4;

%% making a single h time series for ppt

dir_path = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\';
S = dir([dir_path '*.mat']);

ppt_fday = zeros([1,16*86400]);
ppt_h = zeros([1,16*86400]);
ppt_z = ppt_h;
st = 1;
en = 86400;

for i = 1:length(S),
    
    eval(['load ' dir_path S(i).name]);
    ts = data.DATA(:,1);
    fday = data.FDAY;
    
    ppt_fday(st:en) = fday;
    ppt_h(st:en) = ts;
    ts = data.DATA(:,3);
    ppt_z(st:en) = ts;


    st = en+1;
    en = en + 86400;
end;

%% making sinle ts for ppt and ipc

dir_path = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\';
S = dir([dir_path '*.mat']);

ipc_fday = zeros([1,16*86400]);
ipc_h = zeros([1,16*86400]);
ipc_z = ppt_h;
st = 1;
en = 86400;

for i = 1:length(S),
    
    eval(['load ' dir_path S(i).name]);
    ts = data.DATA(:,1);
    fday = data.FDAY;
    
    ipc_fday(st:en) = fday;
    ipc_h(st:en) = ts;
    ts = data.DATA(:,3);
    ipc_z(st:en) = ts;


    st = en+1;
    en = en + 86400;
end;


%% new analysis 


%load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_ppt_z_h_20100218_20100301;
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_ppt_z_h_20100218_20100305;

min_fday = datenum(2010,2,27,0,0,0);
max_fday = datenum(2010,2,27,23,59,59);

L = ipc_fday >= min_fday & ipc_fday <= max_fday;
ipc_fday(~L) = [];
ipc_h(~L) = [];     
ipc_z(~L) = [];     
ppt_fday(~L) = [];  
ppt_h(~L) = [];     
ppt_z(~L) = [];


%b1 = min(ipc_fday):(1/24):max(ipc_fday);
b1 = min(ipc_fday):(2/24):max(ipc_fday);
sp=spline(b1,ipc_h(1:100:end)/spline(b1,eye(length(b1)),ipc_fday(1:100:end)));
ipc_h_f=ipc_h-ppval(ipc_fday,sp);
clear ipc_h;

sp=spline(b1,ipc_z(1:100:end)/spline(b1,eye(length(b1)),ipc_fday(1:100:end)));
ipc_z_f=ipc_z-ppval(ipc_fday,sp);
clear ipc_z;


sp=spline(b1,ppt_h(1:100:end)/spline(b1,eye(length(b1)),ppt_fday(1:100:end)));
ppt_h_f=ppt_h-ppval(ppt_fday,sp);

clear ppt_h;


sp=spline(b1,ppt_z(1:100:end)/spline(b1,eye(length(b1)),ipc_fday(1:100:end)));
ppt_z_f=ppt_z-ppval(ipc_fday,sp);
clear ppt_z;

%clear ipc_z ipc_h ppt_z ppt_h;
%Sc = 1:5:100;
%Sc = 2:5:200;
Sc = 10:20:2000;
%waven = 'morl';
waven = 'cgau4';
perd = 1./scal2frq(Sc,waven,1);


pptzw = cwt(ppt_z_f,Sc,waven);%Complex wavelet transfor
clear ppt_z_f;
ipczw = cwt(ipc_z_f,Sc,waven);%Complex wavelet transform
clear ipc_z_f;
ipchw = cwt(ipc_h_f,Sc,waven);%Complex wavelet transform
clear ipc_h_f;
ppthw = cwt(ppt_h_f,Sc,waven);%Complex wavelet transformf
clear ppt_h_f;

Hxy = ipchw .* conj(ppthw);
L = abs(Hxy) > 100;
Hxy(L) = 100;
ss = sqrt(abs(Hxy));
cc = max(max(ss)) - abs(ss);
weight = (cc./max(max(cc)));
clear cc ss L ipchw ppthw ppt_fday;
% 
%% plot
f = figure;

set(f,'Position',[ 34         562        1876         489]);

b(1) = min(ipc_fday);
b(2) = max(ipc_fday);

% b(1) = datenum(2010,2,27,0,0,0);
% b(2) = datenum(2010,2,28,0,0,0);

ax(1)=subplot(311);
imagesc(ipc_fday(1:10:end),perd/60,abs(ipczw));
set(gca,'FontSize',16);
caxis([0,5]);

% datetick;
axis([b(1) b(2) -inf inf]);
text(b(1),3,'Wavelet amplitude of IPC Z','FontSize', 18,'color','white');
set(gca,'xtick',[]);
colorbar('peer',ax(1),[0.9083 0.7116 0.02452 0.2106],'FontSize',16);

ax(2)=subplot(312)
imagesc(ipc_fday(1:10:end),perd/60,abs(Hxy))
caxis([0,100]);
set(gca,'FontSize',16);
% b=axis;
% datetick;
axis([b(1) b(2) -inf inf]);
text(b(1),3,'Cross wavelet amplitude between IPC H and PPT H','FontSize', 18,'color','white');
set(gca,'xtick',[]);
colorbar('peer',ax(2),[0.9083 0.411 0.02452 0.2106],'FontSize',16);



ax(3)=subplot(313)
imagesc(ipc_fday(1:10:end),perd/60,(abs(ipczw).*weight))
set(gca,'FontSize',16);
caxis([0,5]);
datetick('x','keeplimits');
axis([b(1) b(2) -inf inf]);
text(b(1),3,'Weighted wavelet amplitude of IPC Z ','FontSize', 18,'color','white');
colorbar('peer',ax(3),[0.9083 0.1304 0.02452 0.2106],'FontSize',16);
datetick('x','keeplimits');
xlabel('Date in 2010 mm/dd');

ylabel('Period in minutes');


%% hua processing and plotting al the H components
hua = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\HuaNCAYO_Feb272010.txt');
% DD MM YYYY  HH MM      D       H       Z       I       F
hua_fday = datenum(hua(:,3),hua(:,2),hua(:,1),hua(:,4),hua(:,5),0);


b1 = min(hua_fday):(0.5/24):max(hua_fday);
sp=spline(b1,hua(:,8)'/spline(b1,eye(length(b1)),hua_fday'));
hua_z_f=hua(:,8)-ppval(hua_fday,sp);


load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\ppt20100227vsec.sec.mat;
ppt_z = data.DATA(:,3);
fday = data.FDAY;
ppt_h = data.DATA(:,1);
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\ipc20100227vsec.sec.mat;
ipc_z = data.DATA(:,3);
ipc_h = data.DATA(:,1);
b1 = min(fday):(0.5/24):max(fday);

sp=spline(b1,ipc_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_z_f = ipc_z-ppval(fday,sp);

sp=spline(b1,ppt_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ppt_z_f = ppt_z-ppval(fday,sp);

sp=spline(b1,ipc_h(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_h_f = ipc_h-ppval(fday,sp);
sp=spline(b1,ppt_h(:)'/spline(b1,eye(length(b1)),fday(:)'));
ppt_h_f = ppt_h-ppval(fday,sp);


b(1) = datenum(2010,2,27,6,0,0);
b(2) = datenum(2010,2,27,19,0,0);

%% plot
f=figure;
set(f,'Position',[1          35        1920        1092]);
plot(hua_fday,hua_z_f,'Color',[0 0.749 0.749],'LineWidth', 2);
hold on;
plot(fday(1:10:end),ppt_z_f(1:10:end),'color', [0.6 0.6 0.5],'LineWidth', 2);
plot(fday(1:10:end),ipc_z_f(1:10:end),'r','LineWidth', 2);
datetick('x','keeplimits');
axis([b(1) b(2) -inf inf]);
set(gca,'FontSize',24);
xlabel('Time in hours (UTC)');
ylabel('Z (nT)');
datetick('x','keeplimits');
legend('HUA','PPT','IPM');
kk = datenum(2010,2,27,6,34,0);
%tide_gauge_arrival_fday = 734196.4986110000;%wrong the estimate is 11:50
%UTC
a = axis;
tide_gauge_arrival_fday = datenum(2010,2,27,11,50,0);
fday_eqrthquake = datenum(2010,2,27,06,34,14 );
fday_ttt = datenum(2010,2,27,11,39,0);
 axis([a(1) a(2) -1.5 1.5]);
line([fday_eqrthquake fday_eqrthquake], [-1.5 1.5], 'color',[0.6 0.6 0.6], 'LineWidth', 2);
line([fday_ttt fday_ttt], [-1.5 1.5], 'color',[0.6 0.6 0.6], 'LineWidth', 2);
 text(fday_eqrthquake,-1.2,'Earthquake','FontSize',24);
 text(fday_ttt,-1.2,'Predicted arrival time at IPM','FontSize',24);
 
 % draw box
 
 box_limits = [734196.4547845472   734196.5624171867     -1.       1.];
 line(box_limits([1 1]), box_limits([3 4]),'color',[0.6 0.6 0.6], 'LineWidth', 2);
 line(box_limits([2 2]), box_limits([3 4]),'color',[0.6 0.6 0.6], 'LineWidth', 2);
 line(box_limits([1 2]), box_limits([3 3]),'color',[0.6 0.6 0.6], 'LineWidth', 2);
 line(box_limits([1 2]), box_limits([4 4]),'color',[0.6 0.6 0.6], 'LineWidth', 2);
 
 
%% lOCATION MAP
f=figure;
set(f,'Position',[360        -136        1269         834]);
1
a = imread('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\figures\location_map_150W0N60W60S.PNG');
%a = imread('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\figures\location_map_150W5S60W70S.PNG');
%a = imread('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\figures\location_map_150W0N60W80S.PNG');
a(:,:,1) = flipud(squeeze(a(:,:,1)));
a(:,:,2) = flipud(squeeze(a(:,:,2)));
a(:,:,3) = flipud(squeeze(a(:,:,3)));
%grid2image(a,[18.2,-5,-150]);
grid2image(a,[36.4,-5,-150]);
%set(gca,'Ydir','normal');
set(gca,'FontSize',30);
hold on;
plot(-109.41,-27.17,'b.','MarkerSize',50,'color','blue'); % Easter Islands
plot(-149.57 ,-17.56 ,'b.','MarkerSize',50,'color','blue'); % Tahiti
plot(284.67-360,90-102.05,'b.','MarkerSize',50,'color','blue'); % Huancayo
text(-109.41,-27.17+3,'IPM','FontSize',48,'color','blue'); % Easter Islands
text(-109.41+3,-27.17+3,'Easter Island','FontSize',24,'color','blue'); % Easter Islands
text(-149.57+3,-17.56 ,'PPT','FontSize',48,'color','blue'); % Tahiti
text(284.67-360+3,90-102.05,'HUA','FontSize',48,'color','blue'); % Huancayo
plot(-72.719,-35.846,'w.','MarkerSize',50,'color','white');%epicenter
text(-72.719,-35.846,'Epicenter','FontSize',48,'color','blue');%epicenter
% plot(-125.006,-8.489,'k.','MarkerSize',30,'color','blue'); % DART 51406
% text(-125.006,-8.489,'51406','FontSize',16,'color','blue');% DART
% 17.975 S 86.392 W 32412 

% plot(-86.392,-17.975,'k.','MarkerSize',30,'color','blue'); % DART 32412
% text(-86.392,-17.975,'32412','FontSize',16,'color','blue');% DART 32412
% [rng, az] = distance(-27.17,-109.41,-17.975,-86.392) distance bewteen ipc
% and dart 32412
%[rng, az] = distance(-27.17,-109.41,-35.846,-72.719) dist ipc epicenter
%[rng, az] = distance(-27.17,-109.41,-35.846,-72.719) dist ipc epicenter
%[rng, az] = distance(-27.17,-109.41,-8.489,-125.006) distance bewteen ipc
% and dart 51406
%% plot predicted eta and bz

load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\figures\predicted_eta_data fday_s ipc_z_f_s;
h = 2980; % average water column
Bz = -19017; % Fz from WMM
%h = 3000;

g=9.78; %gravity
c=sqrt(g*h);%phase speed of non-dispersive waves
i_mu = 4*pi*1e-7;%permeabilty 
omega = 4.0;% electrical conductivity sea water (WOA T & S analysis)
K = 1/(i_mu*omega); % magnetic diffusion coefficient
cd = 2*K/h; % lateral diffusion speed
cs = c + sqrt(-1)*cd;% Tylers' complex scaling speed

eta = ipc_z_f_s .* h *cs ./ (Bz *c);
eta_deep = ipc_z_f_s .* h ./ Bz;


plot(fday_s,ipc_z_f_s);
hold on;
plot(fday_s,eta*4,'g','LineWidth',2)
set(gca,'FontSize',16);
datetick;
legend('Observed b_z', 'Predicted \eta')
line([min(fday_s)- 5/(24*60) min(fday_s)- 5/(24*60)], [, -0.5,0.5],'LineWidth', 2)
text(min(fday_s)- 4/(24*60),0,'1 nT','FontSize',16,'color','blue');
line([max(fday_s)+ 5/(24*60) max(fday_s)+ 5/(24*60)], [, -0.25,0.75],'LineWidth', 2,'color','g')
text(max(fday_s)+ 6/(24*60),0.25,'0.25 cm','FontSize',16,'color','green');
set(gca,'ytick',[]);
set(gca,'XMinorTick','on');

%% close up wavelet


load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\figures\predicted_eta_data fday_s
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\ipc20100227vsec.sec.mat;
ipc_z = data.DATA(:,3);
fday = data.FDAY;


b1 = min(fday):(1/24):max(fday);
sp=spline(b1,ipc_z(1:100:end)'/spline(b1,eye(length(b1)),fday(1:100:end)'));
ipc_z_f=ipc_z-ppval(fday,sp);



%waven = 'morl';
waven = 'cgau4';
Sc = 2:10:120;
perd = 1./scal2frq(Sc,waven,5);
ipczw = cwt(ipc_z_f(1:5:end),Sc,waven);

fday_new = fday(1:5:end);

L = fday_new>fday_s(1)-1/24&fday_new<fday_s(end)+1/24;
ipczw(:,~L) = [];
fday_new(~L) = [];

imagesc(fday_new,perd/60,abs(ipczw));
set(gca,'FontSize',16);
datetick('x','keeplimits');
%% making single ts for h, z for ipc and ppt to add data from march 2 to 5



dir_path = 'D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ppt_data\';
S = dir([dir_path '*.mat']);

ppt_fday = zeros([1,16*86400]);
ppt_h = zeros([1,16*86400]);
st = 1;
en = 86400;

for i = 1:length(S),
    
    eval(['load ' dir_path S(i).name]);
    ts = data.DATA(:,1);
    fday = data.FDAY;
    
    ppt_fday(st:en) = fday;
    ppt_h(st:en) = ts;

    st = en+1;
    en = en + 86400;
end;


%% loading abd plotting DART data
data1 = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\DART_32412_Feb27_2010.txt');%g
data2 = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\DART_51406 _feb27.txt');%g
data3 = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\dart51426.txt');%g
%data1 = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\dart51425.txt');
%data1 = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart32411.txt');
%data1 = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart54401.txt');

dart_32412 = data1(:,end);
fday_32412 = datenum(data1(:,1),data1(:,2),data1(:,3),data1(:,4),data1(:,5),data1(:,6));
dart_51406 = data2(:,end);
fday_51406 = datenum(data2(:,1),data2(:,2),data2(:,3),data2(:,4),data2(:,5),data2(:,6));
dart_51426 = data3(:,end);
fday_51426 = datenum(data3(:,1),data3(:,2),data3(:,3),data3(:,4),data3(:,5),data3(:,6));


b1 = min(fday_32412):(0.5/24):max(fday_32412);
sp=spline(b1,dart_32412'/spline(b1,eye(length(b1)),fday_32412'));
dart_32412_f=dart_32412-ppval(fday_32412,sp);

b1 = min(fday_51406):(0.5/24):max(fday_51406);
sp=spline(b1,dart_51406'/spline(b1,eye(length(b1)),fday_51406'));
dart_51406_f=dart_51406-ppval(fday_51406,sp);

b1 = min(fday_51426):(0.5/24):max(fday_51426);
sp=spline(b1,dart_51426'/spline(b1,eye(length(b1)),fday_51426'));
dart_51426_f=dart_51426-ppval(fday_51426,sp);


plot(fday_32412,dart_32412_f,'k');
hold on;
plot(fday_51406,dart_51406_f,'b');
plot(fday_51426,dart_51426_f,'r');
load  D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\offeset_data_for_dart x;
%subplot(411);
plot(fday_32412,dart_32412_f,'k');
a = axis;
hold on;
%subplot(412);
plot(fday_51406 - (x(2,1)-x(1,1)) ,dart_51406_f,'b');
axis(a);
%subplot(413);
plot(fday_51426 - (x(3,1)-x(1,1)) ,dart_51426_f,'r');
%subplot(414);
axis(a);

load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\figures\predicted_eta_data fday_s ipc_z_f_s;
h = 2980; % average water column
Bz = -19017; % Fz from WMM
eta = ipc_z_f_s .* h / Bz;
offset = min(fday_s)-x(1,1);
%subplot(414);
plot(fday_s-offset,eta,'g','LineWidth',1);
axis(a);

% axis([734196.4210349462  734196.7758736559    -0.1424566474     0.2763294798]);
% 
% 
% for i = 1:100,
%     
%     x = ginput(1);
%     offset = min(fday_s)-x(1,1);
%     h=plot(fday_s-offset,eta,'g','LineWidth',1);
%     pause;
%     delete(h(1));
% end;
%     
% 
majs

%% map

worldmap([-45,0],[-180,-60]);
load coast;
geoshow(lat, long,'color','k');
plotm(-27.17,-109.41,'k.','MarkerSize',20,'color','green'); % Easter Islands
plotm(-17.56,-149.57 ,'k.','MarkerSize',20,'color','blue'); % Tahiti
% plotm(284.67-360,90-102.05,'k.','MarkerSize',30,'color','blue'); % Huancayo
textm(-27.17,-109.41+3,'IPC','FontSize',16,'color','green'); % Easter Islands
textm(-17.56 ,-149.57+3,'PPT','FontSize',16,'color','blue'); % Tahiti
% text(284.67-360+3,90-102.05,'HUA','FontSize',16,'color','blue'); % Huancayo
plotm(-35.846,-72.719,'k.','MarkerSize',20,'color','blue');%epicenter
textm(-35.846,-72.719+2,'Epicenter','FontSize',16,'color','blue');%epicenter

plotm(-17.975,-86.392,'k>','MarkerSize',10,'color','black'); %dart 32412
textm(-17.975,-86.392+3,'32412','FontSize',16,'color','black'); % dart 32412
plotm(-8.489 ,-125.006 ,'b>','MarkerSize',10,'color','blue'); % 51406 
textm(-8.489 ,-125.006+3,'51406 ','FontSize',16,'color','blue'); % 51406 
plotm(-22.993 ,-168.098 ,'r>','MarkerSize',10,'color','red');%51426 
textm(-22.993,-168.098 +3,'51426 ','FontSize',16,'color','red');%51426 


%% tide guage data of easter Island Processing and plotting with predicted
%% eta
tide = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\Tidal_Gauge_data_easter.txt');
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\tide_gauge_plot_limits x;

tide_fday = floor(tide(:,1)*10000)/10000;
[a,ia,ib] = unique(tide_fday);
tide_new = tide(ia,:);
tide_new(5:9:end,:) = [];
b1 = min(tide_new(:,1)):(2/24):max(tide_new(:,1));
sp=spline(b1,tide_new(:,2)'/spline(b1,eye(length(b1)),tide_new(:,1)'));
tide_f=tide_new(:,2)-ppval(tide_new(:,1),sp);


load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\ipc20100227vsec.sec.mat;
ipc_z = data.DATA(:,3);
fday = data.FDAY;
b1 = min(fday):(1/24):max(fday);
sp=spline(b1,ipc_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_z_f = ipc_z-ppval(fday,sp);

b1 = min(fday):(1/24):max(fday);
sp=spline(b1,ipc_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_z_f = ipc_z-ppval(fday,sp);
L = fday> x(1,1) & fday < x(2,1);
fday_s = fday(L);
ipc_z_f_s = ipc_z_f(L);
h = 2980; % average water column
Bz = -19017; % Fz from WMM
eta = ipc_z_f_s .* h / Bz;

f=figure;
set(f,'Position',[ 206         246        1040         434]);
plot(tide_new(:,1)+datenum(2009,12,31),tide_f,'b','LineWidth',2);
set(gca,'FontSize',16);
hold;
plot(fday_s,eta,'g','LineWidth',2);
a=axis;

axis([x(1,1) x(2,1) a(3) a(4)]);
datetick('x','keeplimits');
ylabel('Sea level variations (m)');
xlabel('Time in Hours (UTC)');
legend('Tide Gauge','Magnetometer');

%% close up wavelet new

load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\figures\predicted_eta_data fday_s
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\ipc20100227vsec.sec.mat;
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\tide_gauge_plot_limits x;

ipc_z = data.DATA(:,3);
fday = data.FDAY;


b1 = min(fday):(1/24):max(fday);
sp=spline(b1,ipc_z(1:100:end)'/spline(b1,eye(length(b1)),fday(1:100:end)'));
ipc_z_f=ipc_z-ppval(fday,sp);



waven = 'morl';
%waven = 'cgau4';
Sc = 2:10:200;
perd = 1./scal2frq(Sc,waven,5);
ipczw = cwt(ipc_z_f(1:5:end),Sc,waven);

fday_new = fday(1:5:end);

L = fday_new> x(1,1) & fday_new < x(2,1);


ipczw(:,~L) = [];
fday_new(~L) = [];

imagesc(fday_new,perd/60,ipczw./(max(max(ipczw))));
load  d:\manoj\tsunami\COLORTABLE.mat;
colormap(map);
set(gca,'FontSize',32);
datetick('x','keeplimits');
colorbar('peer',gca,[0.9083 0.7116 0.02452 0.2106],'FontSize',32);


%% processing tahiti tide data
tide = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\Tahiti_tide.txt');
tide_fday = floor(tide(:,1)*10000)/10000;
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\Tahiti_tide_data_sel x;

[a,ia,ib] = unique(tide_fday);
tide_new = tide(ia,:);
tide_new(5:9:end,:) = [];

L = tide_new(:,1)> x(1,1) & tide_new(:,1)< x(2,1);
tide_new(~L,:) = [];

b1 = min(tide_new(:,1)):(2/24):max(tide_new(:,1));
sp=spline(b1,tide_new(:,2)'/spline(b1,eye(length(b1)),tide_new(:,1)'));
tide_f=tide_new(:,2)-ppval(tide_new(:,1),sp);
plot(tide_new(:,1)+datenum(2009,12,31),tide_f,'b','LineWidth',2);


tide_e = load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\Tidal_Gauge_data_easter.txt');

tide_fday_e = floor(tide_e(:,1)*10000)/10000;
[a,ia,ib] = unique(tide_fday_e);
tide_new_e = tide_e(ia,:);
tide_new_e(5:9:end,:) = [];

L = tide_new_e(:,1)> x(1,1) & tide_new_e(:,1)< x(2,1);
tide_new_e(~L,:) = [];


b1 = min(tide_new_e(:,1)):(2/24):max(tide_new_e(:,1));
sp=spline(b1,tide_new_e(:,2)'/spline(b1,eye(length(b1)),tide_new_e(:,1)'));
tide_e_f=tide_new_e(:,2)-ppval(tide_new_e(:,1),sp);
subplot(212);
plot(tide_new(:,1)+datenum(2009,12,31),tide_f,'b','LineWidth',2);
hold on;
plot(tide_new_e(:,1)+datenum(2009,12,31),tide_e_f,'r','LineWidth',1);
subplot(211)
a=axis
subplot(212)
axis([a(1) a(2) -inf inf])
set(gca,'FontSize',16);
subplot(211)
set(gca,'FontSize',16);


%% Compare the observed Z with the Seismic data at RPN

load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\seismic_rpn\rpN_Seismic_dta.mat
load D:\Manoj\tsunami\2010_Chile_Easter_Island_data\ipc_data\ipc20100227vsec.sec.mat;
ipc_z = data.DATA(:,3);
fday = data.FDAY;
b1 = min(fday):(1/24):max(fday);
sp=spline(b1,ipc_z(:)'/spline(b1,eye(length(b1)),fday(:)'));
ipc_z_f = ipc_z-ppval(fday,sp);
plot(fday_rpn,rpn_e(:,3)/1e5,'r');
hold on;
plot(fday,ipc_z_f);
datetick;

% The above comparison was done to verify whther the mag signals are
% produced by tilt of the instruments by tsunami waves. The problem : we
% are comparing acceleration , where as tilt is a displacement ? 
% Paper by Yuan and others (Kind gtoup, GFZ, GRL ) has shown that the max
% acceleration for Ind 2004 was due to a tilt of 4e-6 rad (equivalent to
% 0.2 arc seconds that may produce 0.06 nT). So tilt from tsunamis cannot 
% explain the observed1 nT signals. Period ! Sep 6, 2010.

%%
%plot all ipc components - observed and modelled
a =    [734196 + 8/24   734196 + 15/24];
offset_n = 0/(60*24);
 figure(1);
 set(gcf,'Position',[680   119   560   979]);
subplot(311)
plot(fday,ipc_z_f,'LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(fday_simulation + offset_n , ipc_sim_z, 'r','LineWidth',2)
axis([a(1) a(2) -2 2]);
datetick('x','keeplimits');
legend('Observed Z','Predicted Z');
subplot(312);
plot(fday,ipc_x_f,'LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(fday_simulation + offset_n , ipc_sim_x, 'r','LineWidth',2)
axis([a(1) a(2) -2 2]);
datetick('x','keeplimits');
legend('Observed X','Predicted X');
subplot(313);
plot(fday,ipc_y_f,'LineWidth',2);
set(gca,'FontSize',16);
hold on;
plot(fday_simulation + offset_n , ipc_sim_y, 'r','LineWidth',2)
axis([a(1) a(2) -2 2]);
datetick('x','keeplimits');
legend('Observed Y','Predicted Y');


%% wavelet processing of tide gauge data


load('D:\Manoj\tsunami\2010_Chile_Easter_Island_data\dart_tide\Easter_Island_Tide_Digitized.mat');
ts= x(:,2)-nanmean(x(:,2));
waven = 'cgau4';
perdiod_in_minutes = [1:50];

% Sc = scal2frq(1./(fliplr(perdiod_in_minutes)*60),waven,1); %to get desired scales for a set of frequencies
% Sc = scal2frq(logspace(-3,-1,10),waven,1); %to get desired scales for a set of frequencies
Sc = [100:10:1000];

c = cwt(ts,fliplr(Sc),waven);%Complex wavelet transform
perd = 1./scal2frq(Sc,waven,1);
subplot(212);

%imagesc(fday,fliplr(perd)/60,abs(c));
imagesc(fday,fliplr(perd)/60,abs(c(:,1:10:end)));
datetick;
caxis([0,12]);
colorbar('horiz');

